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Analytical Hartree-Fock gradients with respect to the cell parameter have been implemented in 
, the electronic structure code CRYSTAL, for the case of three-dimensional periodicity. The code is 

D ' based on Gaussian type orbitals, and the summation of the Coulomb energy is performed with the 

Ewald method. It is shown that a high accuracy of the cell gradient can be achieved. 
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C/3 ; I. INTRODUCTION 

Electronic structure codes are nowadays widely used 
by theoreticians and also by a growing number of exper- 
y \ \ imentalists. One of the targets is the calculation of the 
• total energy and the structural optimization of large sys- 
^ ' tems. It is well known that this is greatly facilitated by 
I the availability of analytical gradients, and thus the cod- 
^ ing of gradients has become an important part of code 
^ , development. 

^ ' In the molecular case, the geometrical parameters that 
I— 'I must be optimized, are the nuclear coordinates. This 
field was pioneered by Pulay-'^; however it should be men- 
tioned that the theory had already been derived earlier 
independently^. Meanwhile, numerous review articles 
\^ ' about analytical derivatives have appeared'^"^. 
\^ , In the case of periodic systems, the cell dimensions 
CSj ■ are an additional set of parameters. Nowadays, nearly 
all solid-state codes compute the total energy with den- 
sity functional methods. However, because of the suc- 
cess of hybrid functionals, which use an admixture of 
, Fock exchange, there is a growing interest in codes which 

a compute the Fock exchange. CRYSTAL^*^ was originally 
I , a code for pure Hartree-Fock calculations. In the past 
'"^ decade, density functional calculations for all types of 
^ \ functionals have become possible as well with this code. 
O • The code is based on Gaussian type orbitals, and 
most of the contributions to the total energy rely on 
the Hartree-Fock formulation. Therefore, the Hartree- 
Fock gradients with respect to the cell parameter were 
the first step to make gradients with respect to the cell 
' parameter available. In this article, we will report on the 
implementation of gradients at the Hartree-Fock level, 
with respect to the cell parameter, for systems periodic 
in three dimensions. 

The first gradients with respect to the cell parameter, 
at the Hartree-Fock level, were for systems periodic in 
one dimension^^. Meanwhile, various groups have im- 
plemented these gradients in one dimension^^'^^ or in 
two dimensions'^. For the general case, a strategy to 
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compute cell parameter derivatives (and thus the stress 
tensor) was suggested with point charges^^, and an al- 
gorithm for structural optimization, based on redundant 
internal coordinates was proposed^^. 

All these codes use a real space approach, where all 
the summations arc performed in direct space. As an 
acceleration tool, the fast multipole method is applied 
in various cases. The cell parameter gradient is then 
obtained essentially by multiplying the contributions to 
the nuclear gradients with the appropriate factors. 

In contrast, the CRYSTAL code uses, in the case 
of three-dimensional periodicity, the Ewald method 
which is a combination of direct and reciprocal lattice 
summations'^. This means, that besides some contri- 
butions which have to be multiplied with trivial factors, 
there will also be additional derivatives, because various 
parameters in the Ewald scheme, and the reciprocal lat- 
tice vectors, depend on the cell parameters. The Hartree- 
Fock gradients with respect to nuclear coordinates were 
implemented earliei"'^*'''', and after the implementation of 
a tool for structural optimization, it was demonstrated 
that an efficient geometry optimization of large systems 
with any periodicity can be performed^". 

In this article, we describe the formalism used in the 
CRYSTAL code for the cell parameter gradients, and 
present results from tests on various systems. The article 
is structured as follows: first, the variables are defined. 
Then the integrals and their derivatives with respect to 
the cell parameters are discussed, and in the following 
section the derivative of the total energy is given. Fi- 
nally, some examples are used as an illustration of the 
formalism. 



II. VARIABLES 

The primitive cell is given by three vectors: ai, 3,2 and 
03 , and the derivatives with respect to the cell parameters 
aij have been coded, a^- are defined in such a way that 
ail = o-ix is the x-component of ai, a\2 = aiy the y- 
component of a\, and a\z = a\z the 2;-component of ai, 
i.e.: 

ai \ / a\x a\y aiz \ / cm 012 013 \ 

02 = a2x a2y a2z = 021 022 023 (1) 
as / \ a3x a-dy azz J \ 0,31 0,32 0,33 J 

A point g of the direct lattice is defined as ^ = riiai -|- 
71202 + n^o.^, with ni,n2,n3 being integer numbers. We 
assume to have N atoms in the unit cell. The position 
of an atom c in a cell at the origin (i.e. <? = 0) is given 
as Ac = /c,iai + f 0,20,2 + f 0,30,3, and then in cell g the 
position will be: 

= (/c,l+«g,l)ai + (/c,2 + ng,2)a2 + (/c,3 + ng,3)a3 

We have used an additional index, i.e. n^,i means fac- 
tor ni of the lattice vector g. The cartesian t component 
(with t being x, y or z) of the vector Ac + g, indicated as 
Ac,t+9t, is thus 
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As all the integrals depend on the position of the nuclei, 
the derivatives of the nuclear coordinates with respect to 
the cell parameters are required: 

= ^ ifc,m + ng,m)6im5jt = {fc,t + "g,i)<5jt (2) 



with the Kronecker symbol Sjt- 

In the following, we will use the notation Ac to indicate 
the position of atom c. However, we also need to define 
basis functions t/)^ centered at a certain nucleus, where 
fi runs over all the basis functions. For example, basis 
functions fi = 1, ...,5 might be centred at atom 1, basis 
functions fi = 6, ...,15 at atom 2 and so on. It is thus 
trivial to assign a certain atom number c to the basis 
function fi: c = c{^,). We could thus use the notation 
but will instead use the simplified notation A^. 
To avoid confusion, Greek indices are used in this case, 

i.e. = ^c(/j) = ^c- 

The basis functions used are real spherical Gaussian 
type functions, and we will use the notation 4'^{r — A^ — 
9)- 

The crystalline orbitals are linear combinations of 
Bloch functions 

*„(f, k)=Y^ a^„(fc)V'/^(r, k) (3) 

which are expanded in terms of real spherical Gaussian 
type functions 

Mr, k)=J2 Mr- K - 9¥^' (4) 



The spin-free density matrix in reciprocal space is de- 
fined as 

P^,(fc) = 2^a^„(fc)<„(fc)e(ef - e„(fc)) (5) 

n 

with the Fermi energy tp and the Heaviside function 
9. In the case of unrestricted Hartree-Fock (UHF)^^, we 
use the notation 

^l{TS)=Y.<n^^)MrS) (6) 



vI/i(f,fc)=^aJ,JA)V^(f,fc) (7) 

for the crystalline orbitals with up and down spin, re- 
spectively. We define the density matrices for spin up 
and spin down as follows: 

PlAk) =Y,al^{k)a:Uk)&{ep - el{k)) (8) 

n 

and 
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PU^) = E <^U{k)<i{k)Q{eF - ei{k)) (9) 

n 

In the following, P^i, refers to the sum Pj^^ + P^^ 
in the UHF case. The density matrices in real space 
P„n„p;, P^n -P^s - are obtained by Fourier transforma- 
tion of the corresponding reciprocal space quantity. 



III. INTEGRALS AND THEIR DERIVATIVES 



A. Nuclear- nuclear repulsion energy 

The Ewald energy of the nuclear repulsion E^^ is ob- 
tained as 

a,b 

with the nuclear charge Zg, and $(r — Ag) being the 
Ewald potential function as defined in reference 22: 

h K ^ ^ 

= -^ + 11 ^nir- ^^«) + E *^(^"'- ^^«) (11) 

K K 

h are the direct lattice vectors, K the reciprocal lat- 
tice vectors. V is the volume, 7 is a screening parameter 
which was optimized^^ to be 7 = (2.8/y^/^)^, in the 
three dimensional case. The prime in the direct lattice 
summation indicates that the summation includes all val- 
ues of the direct lattice vector h, with the exception of 
the case when \r — Aa — h\ vanishes. In this case, the 
term i — =- is omitted from the sum. In the recipro- 

\r-Aa-h\ 

cal lattice series, the prime indicates that all terms with 
K ^0 are included. 

The error function is defined as 

2 P 

erf(p) = — / expi-u^)du (12) 
Vi" Jo 

We consider the Ewald potential as being dependent 
on the variables Ac , h, V, j and K. The derivative with 
respect to the cell parameters thus requires derivatives 
with respect to the centers Ac and the lattice vectors h 
which are similar to the nuclear gradients and have to be 
multiplied with a trivial factor. In addition, the Ewald 
potential depends on the cell parameters through the 
volume V, the screening parameter 7, and the reciprocal 
lattice vectors K. We obtain: 
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Q^NN 

daij 
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5^NN 
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dAe,t 
daij 


+E 
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5^NN 


dht 
daij 


Q^NN 


dV 




dj 






dKt 


dV 
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daij 


dj 


daij 


dKt 


daij 




dV 




dj 


+E 




dKt 


dV 


dttij 


d^ 


daij 


dKt 


daij 



^dE^ 



Kt, with t = 1,2,3, are the components Ki = 
Kx,K2 = Ky,Kz = Kz of K. The derivatives of the 
parameters V, 7, K with respect to Qui j cirG straightfor- 
ward: 



1. Derivative of the volume 

The volume is obtained as the cross product of the cell 
parameters: 

V = ai{a2 X as) 

Thus, the derivatives -§^, are, for example, obtained 
as: 

dV 

= a2yasz - a.2zasv (14) 



dau 
dV 

^ — = a2za3x - a2xazz (15) 
daiy 

dV 

= aayttiz - a^zaiy (16) 



da^ 



2x 



The remaining components can easily be obtained by 
cyclic permutation: 1 — > 2, 2 — > 3, 3 — > 1, or a; — > y, ?/ — > 

Z, Z ^ X. 



2. Derivative of the screening parameter 

The derivative is straightforward: 

dj_ _ dj_dV_ _ 27 dV 
daij dV daij 3^ daij 



(17) 



3. Derivative of the reciprocal lattice vectors 

The reciprocal lattice vectors K can be expressed as 

= ni6i + 71262 + 71363 (18) 

with the primitive vectors bi of the reciprocal lattice 
defined as: 

61 = -y-a2 X as ; 62 = —as x ai ; 63 = — ai x 02 (19) 

We thus need to evaluate the derivatives of the vectors 
bi with respect to the cell parameters. A few examples 
are given below: 



77 a2za3x - ag^aza; (20) 



and thus 



dbi dV bi 



(21) 



or 



or 



da2x da,2x V 



db2 




dbi dV bi 2tt , , 

+ T7{ -«3. (22) 




(23) 



Again, all the other derivatives can be obtained by 
cyclic permutation. The individual derivatives are thus: 



a.b 



K 



^ = ^ E ZaZ, - E 4f-P(-^(^^^ - - hf) + 7^ E -P (-£ + -^KiA, - A) ] ) (25) 



= E ^"^^17 E 77^ + 7^ - oT + - A.,) exp + ii^(A - A„) (26) 



a/^t 2^ " "y^\^(^2)2 ■ V 27 ■ ^ "''Vy 47 

The partial derivative with respect to the nuclear po- 
sitions is just the nuclear gradient which is already avail- 
able in the code: 

1?^ = - E E(^c . - A„ , - .,) f-^-"^l^'l--f-^^') - 

+ y E ^ (-^) H - ^«)) - e^P - A„))) j (27) 

The derivatives with respect to the direct lattice vec- 
tors h are obtained as: 

^ = - E ^^^^ E(-4^ . - A. . - /^O f ^ " "^1^'^^ " - "^'^ M -P(~7(^ - K- Hf)\ 
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B. Overlap integral 



The fundamental integral is the overlap integral. It is 
defined as 

^^0.,- = jMr- A^)<P. {r - X - g)dh (29) 

Its derivative with respect to the cell parameters is 
thus 

^^^.o,^g _ ( ^^^iO^^g dA^,t ^^^oug djA^^t + gt) \ 



^"•ij ^ V ^X,t ^flij d{A^^t + 9t) daij 



where t represents the summation over x,y,z. 
Exploiting translational invariance, we obtain 



^^liOvg _ ^^tiOug ( dA^^t d{A^^t + 9t) 



Y 

t=l "^l^'t 

dS - 

(31) 

dS - 

The derivative g^""^ is identical to the one which is 
required for the calculation of the gradient with respect 
to nuclear positions and thus it only needs to be multi- 
plied with a trivial factor to obtain the derivative with 
respect to the cell parameters. 



C. Kinetic energy integrals 

The evaluation of the kinetic energy integrals leads to 
a combination of overlap integrals: 

^M3.g = / Mr- X) ['l^f) Mr- X - 9)d'r (32) 

When computing the derivatives with respect to the 
cell parameters, we thus obtain in complete analogy to 
equation 31: 

dT - - 9T - - 
daij ~ dA^j ^^^^ 



D. Nuclear attraction integrals 

The nuclear attraction integrals arc defined as 



N - 

IJ-Ovg 



■^ZaJ Mf-A^)^f-Aa)Mr-A.-9)<i'r (34) 



It has been explicitly evaluated in reference 22. The 
derivative with respect to the cell parameters consists 
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now of the following contributions: there are derivatives 
with respect to three centers A^, A^+g and Aa- Similar 
to the nuclear-nuclear repulsion, derivatives with respect 
to V , 7 and K arc required. These derivatives arc sim- 
ilar to the ones appearing when evaluating the nuclear- 
nuclear gradient. 

The derivatives with respect to the direct lattice vec- 
tors can again be obtained by using the nuclear gradients, 
and the rule: 

The part of the derivative coming from the center {Aj^+ 
g) is obtained in the same way, and the derivative with 
respect to Aa in <I>(r — Aa) is obtained from translational 
invariance. All these nuclear gradients simply need to be 
multiplied with the proper factors, for A^, A^, + g and 

Aa. 

We can thus view the nuclear attraction integrals 
^fiOug ^ ^ function of the variables A^, A^, + g, Aa, 
h, V, 7 and K. As a whole, we obtain therefore for the 
derivative: 



fa,i 



da,j - dA^j dA^,, [U,i+ng,i) y ^^^^^ + ^^^^ J 

Translational invariance can be exploited for the fol- 
lowing contribution: 

j Mr- A,)^^£^^^^-^nj^,Mr- X - g)dh = 



E. Multipolar integrals 

In the expression of the total energy, a term with a 
factor 



(36) 



'(Pc; X) = j Pcimrir- Ac)dh (38) 



appears, with Xp being regular solid harmonics^^. 
The charge Pc{r} is defined as: 

Pc{r^ = - E Pug^.oMr-AMr-A.-9) (39) 

The expression which needs to be differentiated has the 
structure 
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T,<i.s = J Mr- AMr- A, - 9)Xr{r- A^)dh (40) 

We can thus also write 

When computing the derivative with respect to the cell 
parameter, this derivative is required for the expression 
^i^Qvg' ^^'^ only dependence on the cell parameter is via 

the nuclear coordinates and +g. Thus, the deriva- 
tive with respect to the cell parameters a^- is obtained 
as 



F. Field integrals 

They are defined as 

Ko^Sc = ^"(^^) J ^Mr- A. - g) ^{?- Ac)-Y.TZ 



pen 

' d^r (43) 

. , Ac-h\ 



with ZY^{Ac) being the spherical gradient operator 
(Ref. 22). The penetration depth pen is a certain thresh- 
old for which the integrals are evaluated exactly^^'^^. 

Similar to the nuclear attraction integrals, this inte- 
gral also requires derivatives with respect to A^, A^ + g 
and Ac, h and the derivatives with respect to V , 7 and 
K because of the Ewald potential. The derivatives with 
respect to A^^ and Ai, + g arc already available, and the 

derivatives with respect to Ac are obtained from transla- 
tional invariance. We obtain 

lnOvgc lnOi^gc n . IfxOvgc / / , \ / IfiOvgc . lnOvgc \ r 

dai^ - ~^A~^^'' + dA^j ^^"'^ + ~ [~dA~ + ~dA~ ) 
+zriAc) y Mr-A,)Mr-A.-g)[j: ' %^ -n.-Y.^.V^T^M'''^ 

dV aij duij ^ dKt daij 



Similar to the nuclear attraction integrals, we can ex- 
ploit translational invariance: 

J A,)Mr- X - g)^^^^^^nj;.dh = 

- / Mr- A,) ^'^'^^ " ^{r- Ac)n^,Ah (45) 
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and 



1 



7^ 



b„{f- A„ - g) 



\r - Ac-h\ 



G. Spheropole 



The spheropole is a correction required because the 
charge distribution is approximated by a model charge 
distribution in the long range^^. 



— V T p - 



(j!)^(f - l^)^!>^(f - X - g)+ 

J Mr' - A^)Mr' - X - 9)Xrir' - A^)d^^ 'SpiA^, 

6p{Ac,r) is a Gaussian representation of a unit point 
multipole^^. 

There are thus dependencies on A^j,, A^ + g, and the 
volume V. We need to compute 

E T.P^S.3^Z J [-Mr-A,)Mr-A.-9)+ 



If^^d^r (47) 



J - A^)Mr' - X - 9)xr{r' - K)^^^ '^T^K^ ^ 



2tt d 



i)f,{f- A^)(j)^{r- A^-g) + 



J Mr' - \)Mr' - X - 9)Xr{r' - A^)dh '6r{A^, r) 



Q dV 
V dan 



(48) 



H. Bielectronic integrals 



These are integrals of the type 



jjQvgrnan+h 



/ 



i>n{r-A^)(j)^{f-A^ - g)(f)rf-Ar - n)Mr'-Aa -n-h) 3 , 



r — r ' 



-d^rd^r' (49) 



The derivative with respect to the cell parameters aij 
is straightforward, as the only dependence on the cell 
parameter is via the position of the c;enters and thus the 
nuclear gradients only need to be multiplied with a factor: 

dB dB ^ dB 



^ IJ,Oi/grfi<Tn-\-h 

dA~ 



OB ^ 



{fT,i + nfi,i) + — (/a,i + nn,i + uj^^^) (51) 
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As in reference 19, we define a Coulomb integral as 

follows: 

pen 

^ /j.OvgrOcTh / ^ fj.Oi/gTnan-\-h 

(52) 

n 

and an exchange integral: 

(53) 

n 

Using this notation, one summation over the lattice 
vectors is already included. 

I. Energy weighted density matrix 

The contributions of the energy weighted density ma- 
trix to the derivative with respect to : 



dAi Jbz 



exp(ifcg) { alj{k)a*^im']{k) + Q)e(e^ - e]{k) - Q) 



+<,{k)a*^,{k){e\{k) + Q)Q{ep - ej(fc) - Q) } d^k 

/ ••• (54) 

JBZ 



require a different prefactor, compared to the case of 
the gradients with respect to nuclear positions: 

oaij j^z Jbz 



IV. TOTAL ENERGY AND GRADIENT 

A. Total energy 

The total energy, in the notation of reference 19, is 
obtained as follows: 

^total ^kinetic ^ j^NN _|_ ^coul — nuc _^ ^coul— el _|_ ^exch— el 

= ^ Z^Z.^iA -A,)+J2 Kg^oT^O.g 

a,b 9,^5^ 



+ 2 E ^'^gtioi ~ Q^iJ.Oug+ E ^ahrQ^ tiOvgrOah ~ EE E ^("(^cl X)-M™g^^^ J 
g,IJ,,v ^ h,T,iT '=0 rn=-l 

_i Y P\^y . . . . - i y - y ^ n (56) 

2 .iL^ ffluO cr/trO fiOiygrOa-h o i—^ ^SmO crhrO fiOiygrOirh V / 



2 .il-^ i^s/iO cr/trO IJ-^vgr^crh 2 .^--^ "S/tO 
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B. Gradient of the total energy 



The gradient with respect to the cell parameter is a 
combination of the nuclear gradients, with appropriate 
factors, and new derivatives with respect to parameters 
such as V, 7 and K and their derivatives with respect to 
the cell parameters. 



^^total 

dan 



da. 



?! II 1/ ^ r' T ^ J*' 



L 



1 n J 



^=0 m=-l 



L I 



p "^^fiOvgrQah \ ^ \ ^ \ ^ \ ^ p 



=0 m=-iS,Tec,(7 



r{f- Ar)Mr- A. - h)Xr{?- A,)d\ M™„.^.^ 



} 

, 1 pT pT ^-^nQj^grOah , 1 y^ pi y^ pi 



-E 



^^ix5vg 



da. 



l^iOvgrOah 
daij 

exp(ifc5) J2 { «L(fcX(^)(4(^) + - eiik) - Q) 



•ij JBZ 



+4n(fc>l(fc)(4(fc) + Qmep - ei{k) - Q) }d^k 



V. EXAMPLES 



In this section, wc give some numerical examples of 
the accuracy of the gradients. First, we consider MgO, 
at a lattice constant close to the equilibrium lattice con- 
stant. In table I, numerical and analytical gradients are 
compared, for various values of the "ITOL" parameters 
controlling the accuracy of the calculation of the inte- 
grals. As was explained in reference 18, certain parame- 
ters (IT0L2, IT0L4 and IT0L5) can introduce an asym- 
metry in the evaluation of the integrals, which results in 
inaccuracies in the gradients. Exactly the same holds for 
gradients with respect to the cell parameter. The ac- 
curacy for the default parameters is about 2x10^^ a.u. 
which should be good enough for practical purposes, and 
by increasing the values of these parameters, the error is 
reduced to 10~^ a.u. 

In table II, various MgO supercells from 1x1x1 (i.e. 
with one magnesium and one oxygen atom in the primi- 
tive cell) up to 5 X 5 X 5 have been considered (i.e. with 
125 magnesium and 125 oxygen atoms in the primitive 
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cell). The results demonstrate the high stability of the 
gradient when larger cells are used. 

In table III, further examples illustrate the accuracy of 
the gradients. For various systems, including magnetic 
ones, the analytical and the numerical gradients agree 
reasonably well. It is again demonstrated that increasing 
the ITOL-parameters leads to more accurate gradients. 
Also, the stability with respect to the supercell size is 
illustrated, for AI2O3. 

In table IV, the total energy and the analytical gradi- 
ent are displayed, around the equilibrium structure. We 
note that in all cases, the gradient changes sign around 
the equilibrium: for example, for MgO, the energy has 
its minimum between 4.18 and 4.19 A, and also the an- 
alytical gradient changes its sign. Similar, in the other 
systems considered, the minimum of the energy and the 
zero of the analytical gradient agree to 0.01 A, at least. 

This is also demonstrated in table V, where the ge- 
ometry has been optimized according to the minimum 
in energy, or the vanishing of the cell gradient (i.e. the 
minimum of the energy was determined, up to an ac- 
curacy of 0.001 A, and similarly, the geometry with the 
smallest value for the gradients was determined, up to 
an accuracy of 0.001 A). It turns out that the two min- 
ima differ at most by 0.004 A, which is probably lower 
than the noise by the other parameters (basis set, choice 
of FIXINDEX parameter^^ and so on). Note that these 
calculations were done with the fractional coordinates of 
the atoms held fixed, by simply varying the cell param- 
eters (an automatic optimizer which optimizes the cell 
dimensions and the nuclear positions simultaneously, us- 
ing analytical gradients, is not yet implemented in the 
CRYSTAL code). 

Finally, in table VI, the CPU times are displayed. The 
calculations were performed on a single CPU of a Com- 
paq ES45, with a clock rate of 1 GHz. It is probably 
best to compare the CPU time for the integrals with the 
time for the gradients, as the code is somewhat similar 
for these two tasks. At present, the CPU time for all 
the gradients (nuclear and cell gradients) is roughly ten 
times the CPU time for the integrals. This ratio is ex- 
pected to be the upper limit as the gradient code is not 
yet fully optimized. However, the calculation of numeri- 
cal gradients scales with the number of parameters to be 
optimized, because at least one more energy point is nec- 
essary for one additional numerical derivative. Thus, if 
there are enough geometrical parameters, the analytical 
gradients should be clearly favorable. 

For the MgO supercells, one can also analyze the CPU 
times for the integrals and the self-consistent field proce- 
dure as a function of the system size. When dividing by 
the number of iterations (which is 14, 14, 15, 15 and 18 
for the cells from size 1x1x1 up to 5x5x5), the CPU 
time per iteration scales roughly with the third power of 
the system size which is to be expected as the diagonal- 
ization scales with this power. The integrals scale with 
a somewhat lower exponent (less than two), due to the 
fact that more and more of the bielectronic integrals of 
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the larger cells are not evaluated exactly, but with the 
help of a multipolar expansion. 

VI. CONCLUSION 

A formalism for the calculation of the analytical gra- 
dient of the Hartrce-Fock energy, with respect to the cell 
parameter, has been presented and implemented in the 
code CRYSTAL, for the case of systems periodic in three 
dimensions. The implementation includes the cases of 
spin-restricted and unrestricted polarization. It has been 
shown that a high accuracy can be achieved. Future de- 
velopments such as a full structural optimization with 
the help of analytical gradients now become feasible. 
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TABLE I. Fee MgO, at a lattice eonstant of 4.25 A. The accuraey of the analytical gradient as a function of the truncation 
parameters ("ITOL" -parameters) is displayed. Basis sets of the size [3s2p] were used for Mg and O. 



ITOL 

6 6 6 6 12 (default) 
8 8 8 8 14 
10 10 10 10 16 
10 10 10 16 16 
10 12 10 16 16 



analytical derivative 
[Eh/ao] 
-0.012737 
-0.012589 
-0.012522 
-0.012496 
-0.012505 



numerical derivative 
[Eh/ao] 
-0.012555 
-0.012533 
-0.012471 
-0.012482 
-0.012503 



TABLE II. Fee MgO, as in table 1. The analytical gradient is computed as a function of the supercell size to demonstrate 

the numerical stability (from 2 to 250 atoms per cell). The default ITOL-parameters (6 6 6 6 12) are used. 

supercell size total energy/MgO unit analytical derivative 

[Eh] [Eh/ao] 

1x1x1 -274.6635207 -0.01273658 

2x2x2 -274.6635204 -0.01273664 

3x3x3 -274.6635205 -0.01273668 

4x4x4 -274.6635204 -0.01273665 

5x5x5 -274.6635204 -0.01273665 



TABLE III. Other examples for a comparison of analytical and numerical gradient, including ferromagnetic (FM) and 
antiferromagnetic (AF) states. If not stated otherwise, the default ITOL parameters are used. The basis sets are in the range 
from [2slp] for H in urea, up to [5s4p2d] for the transition metals. 



system 
AI2O3 
AI2O3 

Urea 

Urea 
NiO, FM 

NiO, FM (ITOL: 10 12 10 16 16) 225 
NiO, AF 225 
NiO, AF (ITOL: 10 12 10 16 16) 225 
KMnFs, FM 221 



space cell parameters component 

group [A] 
167 4.7602, 12.9933 --^ 



167 4.7602, 12.9933 
113 5.565, 4.684 
5.565, 4.684 
4.20 



113 

225 



4.20 
4.20 
4.20 
4.19 



as 



dE 



dE 
aiz 



dE 
a3z 



dE 



dE 



analytical derivative 

[Eh/ao] 

-0.19630 (2x2x2 cell: -.19630) 
-0.06366 (2x2x2 cell: -.06366) 

-0.01501 

-0.02495 

0.00595 

0.00591 

0.01111 

0.01094 

0.01043 



numerical derivative 

[Eh/ao] 
-0.19625 

-0.06361 

-0.01475 

-0.02516 

0.00656 

0.00592 

0.01234 

0.01109 

0.01095 
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TABLE IV. Other examples for a comparison of analytical and numerical gradient. The default ITOL parameters are used. 



system 


cell parameter 


energy 


components 


gradient 




[A] 


[Eh] 






MgO 


4.18 


-274.664192 


d 


8.495x10""' 




4.19 


-274.664222 


8.103x10"® 




4.20 


-274.664209 




-6.735x10""* 


Urea 


5.52, 4.64 


-447.683214 


a a 


7.057x10""*, 1.4379x10"^ 




5.53, 4.63 


-447.683158 




-4.045x10"", 6.1033x10"^ 




5.53, 4.64 


-447.683218 




-7.904x10"^*, 6.649x10"" 




5.53, 4.65 


-447.683176 




-1.1725x10"^, -4.7011x10"^ 




5.54, 4.64 


-447.683166 




-2.2707x10"*, -1.010x10"" 


KMnFg 


4.28 


-2047.643166 


a 

aaix 


4.098x10""* 




4.29 


-2047.643181 


-5.754x10"" 




4.30 


-2047.643141 




-1.5369x10"* 



TABLE V. Optimized structures, using energies or analytical gradients. The default ITOL parameters are used. For each 
compound, the upper line refers to the structure with the lowest energy, and the lower line to the structure with (practically) 



vanishing 


force. The components of the forces 


are as in table III. 




system 


geometry 


energy 


force 




[A] 


[Eh] 


[Eh/ao] 


KMnFg 


4.288 


-2047.643182 


-3.8-10"" 




4.284 


-2047.643179 


1.3-10"® 


Urea 


5.525 ; 4.642 


-447.683224 


-1.2-10"" ; -2.9-10"® 




5.524 ; 4.642 


-447.683224 


2.8-10"® ; 4.8-10"® 


AI2O3 


4.497 ; 12.111 


-1401.048515 


-4.4-10"" ; -1.3-10"" 




4.496; 12.111 


-1401.048515 


2.9-10"" ; -2.2-10"® 



TABLE VI. CPU times for the various calculations. The calculations were performed on a Compaq ES45, using a single 
CPU (1 GHz). The CPU times refer to the part for the integrals (all the integrals were written to disk), the self-consistent 
field (SCF) procedure, and to the calculation of all the gradients (i.e. nuclear gradients and cell gradients). 



system 


number of 




CPU time, in seconds 






symmetry operators 


integrals 


SCF 


gradients 


MgO (1x1x1) 


48 


2 


0.5 


26 


MgO (2x2x2) 


48 


11 


18 


152 


MgO (3x3x3) 


48 


55 


500 


533 


MgO (4x4x4) 


48 


209 


6330 


1662 


MgO (5x5x5) 


48 


670 


57851 


4443 


AI2O3 (1x1x1) 


12 


15 


10 


184 


AI2O3 (2x2x2) 


6 


544 


4681 


3877 


Urea 


8 


29 


103 


257 


NiO, FM 


48 


12 


6 


128 


NiO, AF 


12 


32 


220 


346 


KMnFa 


48 


27 


20 


281 
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